/*
  Copyright 2010-2012 Patrik Jonsson, sunrise@familjenjonsson.org

  This file is part of Sunrise.

  Sunrise is free software; you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation; either version 3 of the License, or
  (at your option) any later version.

  Sunrise is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with Sunrise.  If not, see <http://www.gnu.org/licenses/>.
  
*/

/** This tests that the output from thermal_equilibrium_grain with
    bbdust.fits (no scattering, flat opacity) is consistent with what
    you get from the analytic blackbody_grain_model. */

#include<iostream>
#include "CCfits/CCfits"
#include "counter.h"
#include "polychromatic_grid.h"
#include "optical.h"
#include "scatterer.h"
#include "constants.h"
#include "dust_model.h"
#include "emission.h"
#include "emission-fits.h"
#include "full_sed_emergence.h"
#include "full_sed_grid.h"
#include "emergence-fits.h"
#include "mcrx-units.h"
#include "grid_factory.h"
#include "grid-fits.h"
#include "shoot.h"
#include "fits-utilities.h"
#include "boost/shared_ptr.hpp"
#include "density_generator.h"
#include "preferences.h"
#include "model.h"
#include "wd01_Brent_PAH_grain_model.h"
#include "wd01_grain_model.h"
#include "xfer.h"
#include "ir_grid.h"
#include "equilibrium.h"
#include "vecops.h"
#include "blackbody_grain.h"
#include "misc.h"

using namespace mcrx;
using namespace CCfits;
using namespace std; 
using namespace blitz; 

int main(int, char**)
{
  T_unit_map units;
  units ["length"] = "m";
  units ["mass"] = "kg";
  units ["luminosity"] = "W";
  units ["wavelength"] = "m";
  units ["L_lambda"] = "W/m";

  cout << "This tests that the output from thermal_equilibrium_grain with bbdust.fits\n(no scattering, flat opacity) is consistent with what you get from the\nanalytic blackbody_grain_model.\n\n";

  const int n_threads=1;

  array_1 lambda(logspace(1e-7,1e-3,150));

  Preferences p;
  p.setValue("grain_data_directory", string(word_expand("$HOME/dust_data/crosssections")[0]));
  p.setValue("wd01_parameter_set", string("DL07_MW3.1_60"));
  p.setValue("use_dl07_opacities", true);
  p.setValue("n_threads", n_threads);
  p.setValue("use_grain_temp_lookup", false);

  blackbody bb(5800,1.55e-9*1000000);
  array_2 intensity(1,lambda.size());
  intensity(0,Range::all())=bb.emission(lambda);

  thermal_equilibrium_grain teg(word_expand("$HOME/dust_data/crosssections/bbdust.fits")[0],p);
  cout << "bbdust has cross section " << teg.sigma(0)(0) << endl;

  // create a blackbody grain with *opacity* equal to the bbdust cross section
  blackbody_grain_model<polychromatic_scatterer_policy, mcrx_rng_policy> 
    bbg(teg.sigma(0)(0), units);
  teg.resample(lambda);
  bbg.set_wavelength(lambda);

  array_2 crapsed(2,lambda.size());
  crapsed=0;
  array_2 crapsed0(crapsed(Range(0,0),Range::all()));
  array_2 crapsed1(crapsed(Range(1,1),Range::all()));
  array_1 crapmdust(1);
  crapmdust=1;

  // to get the equivalent amount of emission from the
  // thermal_eq_grain, we supply mdust=1 and dn=1, which means that
  // the opacity we get is equal to the cross section, which is what
  // we specified for the blackbody grain.
  teg.calculate_SED_from_intensity_virtual(intensity, crapmdust,crapmdust,
					   crapsed0,false);
  bbg.calculate_SED_virtual(intensity,crapmdust,crapsed1);

  ofstream of("sed.txt");
  for (int i=0; i<lambda.size(); ++i)
    of << lambda(i) << '\t' << crapsed0(0,i) << '\t' << crapsed1(0,i) << '\n';


  const T_float teg_lint= integrate_quantity(crapsed0(0,Range::all()),
					     lambda,false);
  const T_float bbg_lint= integrate_quantity(crapsed1(0,Range::all()),
					     lambda,false);

  cout << "Low temp:\n";
  cout << "Emitted luminosities: " << teg_lint << '\t' << bbg_lint << endl;
  const T_float dev1=teg_lint/bbg_lint-1;
  cout << "Deviation: " << dev1 << endl;
  cout << "Low temp, max deviation is " << max(crapsed1-crapsed0) << endl;

  crapsed=0;
  teg.calculate_SED_from_intensity_virtual(1e4*intensity, crapmdust,crapmdust,
					   crapsed0,false);
  bbg.calculate_SED_virtual(1e4*intensity,crapmdust,crapsed1);

  const T_float teg_lint2= integrate_quantity(crapsed0(0,Range::all()),
					      lambda,false);
  const T_float bbg_lint2= integrate_quantity(crapsed1(0,Range::all()),
					      lambda,false);
  const T_float dev2=teg_lint/bbg_lint-1;

  cout << "High temp:\n";
  cout << "Emitted luminosities: " << teg_lint2 << '\t' << bbg_lint2 << endl;
  cout << "Deviation: " << dev2 << endl;

  const T_float tolerance=1e-5;
  if((fabs(dev1)>tolerance) || (fabs(dev2)>tolerance)) {
    cerr << "\nFAILED!\n\n";
    return 1;
  }
  else {
    cout << "\nPassed\n\n";
    return 0;
  }
}
